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ABSTRACT 

This paper outlines a number of finite element methods for the analysis of 
nonlinear problems in rubber elasticity. Several different finite element schemes 
are discussed. These include the augmented Lagrangian method, continuation or 
incremental loading methods, and associated Riks-type methods which have the capa- 
bility of incorporating limit point behavior and bifurcations. Algorithms for the 
analysis of limit point behavior and bifurcations are described and the results of 
several numerical experiments are presented. In addition, a brief survey of some 
recent work on modelling contact and friction in elasticity problems is given. 
These results pertain to the use of new nonlocal and nonlinear friction laws. 
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PRELIMINARIES AND NOTATIONS 


KINEMATICS 


The usual notation in finite elasticity is employed: U is the displacement 

vector, X is the position of a particle in the current configuration whose posi 
tion was X in the reference configuration, F is the deformation gradient. 
These quantities are illustrated in the figure~below. 


' * • . • , 
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(x'x’X) uj 


DISPLACEMENT VECTOR 
X - X 


DEFORMATION GRADIENT 


V(X + u) ; F - = Sf*. {X A 
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HOMOGENEOUS, ISOTROPIC 
HYPERELASTIC MATERIAL 

We begin with a study of equilibrium problems in finite elasticity. It is 
assumed that the materials involved are hyperelastic, isotropic, and homogeneous. 
Therefore, they are characterized by a strain energy function W which is given 
as a function of invariants of the deformation tensor Q . These problems are 
complicated by the fact that nonconvex constraints must~be enforced. For compres- 
sible materials the constraint manifests itself in the condition that J(u) = det 
F > 0 while for incompressible materials J(u) =1. 


W = W (I-,,I 2 ,J) = STRAIN ENERGY PER UNIT VOLUME 
^ = TRACE C 

l 2 = \ (trace C) 2 - Jr trace c 2 
J = Y^3 = det F C = F T F 

COMPRESSIBLE MATERIALS 

J(u) = det F > 0 

INCOMPRESSIBLE MATERIALS 

J(u) = i 
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SOME CONDITIONS ON THE FORM OF W 

Some conditions on the form of the energy W are presented below. In addi- 
tion to the fact that the energy must be form invariant under changes of the 
spatial frame of reference, other conditions must be enforced if one expects a 
well-behaved solution. Three of these are listed below: 


1. COERCIVITY : 

w(u) = w(i 1 (u),i 2 (u),j(u)) 

W(u) -* +°° as j | Vu 1 1 -+ +°o 


2. SINGULAR BEHAVIOR: 


W(u) -*■ +°° 


aw 

3F 


-> +00 


as det F -+ 0 


3. QUASICONVEXITY: 


a w(u) 
a f . a . 


A-A.y y > 0 
J a B — 


V X , u w 
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INCOMPRESSIBLE AND NEARLY INCOMPRESSIBLE MATERIALS 

« 

Two basic classes of methods are employed here to handle compressible and 
nearly incompressible materials: Lagrange-multiplier methods (mixed methods), in 
which the incompressibility constraint h(j) = 0 is accounted for using Lagrange 
multipliers and penalty methods in which the total potential energy func- 
tional II is penalized by the addition of a positive semi-definite, generally 
convex penalty functional. 


INCOMPRESSIBILITY CONSTRAINT 

h ( J ) = 0 

h { J ) = J - 1 , J 2 - 1 , (J - l) 2 , -InJ , etc. 
W = W ( I n , 1 2 ) - p h ( J ) 

p = Lagrange Mult. -Hydrostatic Pressure 


PENALTY TERMS 


n(u) = 


W(u)dX - f(u,v) 


0 . 


+ e 


-1 


'ft 


g(J(u))dX 


g (J ) > o g ( J ) = 0 J = 1 
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EQUILIBRIUM CONDITIONS 


The finite elements employed are based on various alternative statements of 
equilibrium conditions for elastic bodies: 

1) ENERGY FORMULATION 

n(u) = TOTAL POTENTIAL ENERGY 

W(u)dX + f ( u , V ) 
fl - ” " 

f(u,v) = f f(u)*vdX + f t(u)*vds 
J o' ~ j r 2 ~ ~ ' 

JT(u)<n{v)forall v in 

K = { v : | f W ( v )dX | <“ ; 

J 0 “ 

v = 0 on T-j ; det 7v > 0 or 

2) AUGUMENTED LAGRANGE 

L ( u , F , A ) = n(u) + I I |7u - F I 2 dX 
^ J n = 

+ f X : (7u - F ) dX 
Jfi- ' 

L : V x (G: det G = 1 } x { ■ ■ }* 


The augmented Lagrange method combines features of both penalty and 
Lagrange-multiplier schemes. The incompressibility constraint can be satisfied 
a priori in a straightforward manner for Mooney-Rivlin materials, with the 
result that the method is extremely fast and efficient when used in conjunction 
with, e.g., Uzawa's method (ref 1). (See equation.) 
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FINITE ELEMENT METHODS 


SOME ELEMENT FAMILIES 

The figure below illustrates some of the standard finite el em ent methods 
employed. Unfortunately, not all of these methods are numerically stable since some 
may not satisfy the generalized LBB condition of Oden and LeTallec (refs. 2 and 3). 
This condition is given In the equation on the following page, where Q is the space 
of Lagrange multipliers, is the space of finite element approximations of dis- 

placement, p is the hydrostatic pressure and ||Vv h || 0 ,p is an appropriate energy 
norm on v ■ . The existence of a Bh > 0 independent~of mesh size h is necessary for 









The numerical stability of mixed- and penalty finite element 
methods depends upon this condition, and particularly on the behavior 
of the parameter with the mesh size h. 
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STABILITY RESULTS 

Oden and Jacquotte (ref, 4) have recently completed a mathematical analysis of 
the LBB condition for incompressible viscous flows. Some of these results appear to 
be directly applicable to the finite elasticity problem. The following table summar 
izes the behavior of certain finite element methods for these constrained problems. 
The results generally fall into three categories: the stability parameter 3h is 

independent of h and the method is stable, 8ft is dependent on h and the method 
is unstable having spurious pressure modes, or the element is "locked", meaning that 
the penalty parameter e depends on the mesh size and that the displacements 
approach 0 as e tends to zero for a fixed mesh size h . 


VELOCITY APPROX. 


QUADRATURE RULE 
(PRESSURE APPROX. Q ) 


CONVERGENCE 

RATE 












RATE OF CONVERGENCE 


UNSTABLE PRESSURE 


UNSTABLE PRESSURE 


SUBOPTIMAL (0(h)) IN VELOCITY 
ERROR IN ENERGY NORM 


LOCKS FOR SMALL t , t MUST 
BE TAKEN AS DEPENDENT ON h 


RATE OF CONVERGENCE 


UNSTABLE PRESSURE 


OPTIMAL : 


SUBOPTIMAL 10(h ) ) IN VELOCITY 
ERROR IN ENERGY NORM 


UNSTABLE PRESSURE' 

















VELOCITY APPROX. 




0(1) OPTIMAL 


0<fc) UNSTABLE PRESSURE 


x 


0 ( 1 ) 


SUBOPT I MAL (0(h) IN VELOCITY 
ERROR IN ENERGY NORM 












ALGORITHMS 


Several different algorithms are employed for the analysis of finite elasti- 
city problems. These include the augumented Lagrange/Uzawa methods, continuation 
(incremental loading) methods, and homotopy methods. In the present paper, augmented 
Lagrange methods are discussed briefly, but the focus is on the continuation- type 
methods of, for example, Riks (ref. 5), Keller (ref. 6), Crisfield (ref. 7), and 
Padovan (ref. 8). 

1. AUGUMENTED LAGRANGE/UZAWA 

1. EXTREMELY SIMPLE & FAST 

2. FOLLOWS STABLE BRANCHES 

3. NOT EASILY ADAPTED TO GENERAL MAT'LS 

2. CONTINUATION (INCREMENTAL LD'G) 

1. Riks (ref. 5), Wempner (ref. 9) 

2. Keller (ref. 6), Rheinbolt (ref. 10) 

3. Crisfield (ref. 7), Padovan (ref. 8) 

3. HOMOTOPY METHODS 

f(x,p) = o 

3f (x(s ) ,p(s ) ) . 9 f ( x ( S ) , p ( s ) ) . 

2 X + 2 P = 0 

3x ~ 3 p 

N ( x ( s ) , p ( s ) ) = 0 

+ ODE SOLVER 
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AUGUMENTED LAGRANGE METHOD 
The augumented Lagrange is a super-fast method. 


AUGUMENTED LAGRANGE c=C> BLOCK RELAXATION 

1. x n+1 = x n - p ( Vu n - F n ) 


2 . 


L(u£ * Ek+1 5 ~ n) - L( ~k • 9 ’ *") 
| SOLVED EXPLICITLY, N + 2,3| 


3 * 3 v L ^ ~k+l ’ ~k + l * ~ ^ * X 0 


CONTINUATION METHODS 


f ( x ( s ) , p ( s ) ) = 0 

x ( s ) . x(s) + p 2 (s) = 1 £= N ( x , p ) = 0 
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NUMERICAL EXAMPLES 


Consider the inflation of a thick rubber spherical shell subjected to external 
loading. The shell is pressurized until it snaps through. Pressure continues to 
increase until the shell inflates. The material is assumed to be an incompressible 
Mooney-Rivlin material. The numerically stable Q2/P1 element is employed. Good 
results are obtained using the Riks-Crisf ield method (refs. 5 and 7) with Newton- 
Raphson correction. Geometry of the shell is shown below and this is followed by 
several figures which illustrate numerical results. It is noted that the stiffness 
of the shell is strongly dependent on the material properties. In particular, for 
a fixed value of the Mooney-Rivlin constant C^, if C£ is chosen sufficiently small, 
a limit point type behavior occurs in the inflated shell. This represents the phe- 
nomenon of a large decrease in pressure in an inflated balloon with a large increase 
in strain at the crown. This limit point behavior disappears for larger values of 
C2- (See ref. 11.) 
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frb/in^ 

P= 0. 

P= 1.449 
P* 0.52? 

P= 1.591 
Ps14.437 
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CALCULATION OF BIFURCATION PATHS 


The Riks method (ref. 5) makes it possible to calculate bifurcation paths as 
well. A second example involves the buckling of a thick rubber plate under end thrust. 
This problem was analyzed by Sawyers and Rivlin (ref. 12), and provides a good example 
to test the bifurcation capabilities of the code. For low initial aspect ratios, a 
barreling mode of bifurcation is obtained, whereas for thinner slabs a flexural mode 
is obtained. The following figures show computed numerical results for these ele- 
ments. Bifurcation branches were calculated for this problem as well. The fol- 
lowing figures show the results of these calculations. After checking the energy 
of the system on each branch, the branch of lower energy is then calculated. In 
the first figure, branch one has lowest energy and branch two indicates a second- 
ary bifurcation. The lower branches in each of the subsequent figures also rep- 
resent lower energy equilibrium states. (See tef. 13.) 


PLATE UNDER END THRUST 
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first three buckling modes eor a 


PLATE UNDER END THRUST 
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ORIGINAL PAGE IS 
» POOR QUALITY 
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COMPUTATIONAL PROBLEMS DURING HIGH COMPRESSION 


The augumented Lagrange method failed for problems with high compression as 
indicated in the figure below. An oscillatory mode was obtained at compressions of 
30%. Here the Q1/P0 element was employed. In view of this, the continuation tech- 
niques with Q2/P1 elements were attempted in the following examples. 


30% COMPRESSION 
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FIXED END CYLINDER 


Area shown in 
deformation plots 



- M 


- N 












COMPRESSION OF A ROBBER BLOCK 

Here a rubber block of revolution is one-quarter of the cross section of a 
block of revolution is shown. The block is subject to compression by uniformed dis- 
placement of the edges. A strong singularity is developed in the free edge leading 
to a cusp as indicated in the calculated deformed configurations. Ultimately, contact 
of these surrounding regions occurs. This is a very difficult class of elastostatic 
problems and contact conditions must be incorporated in the analysis procedure to 
handle these problems. 


PLANE STRAIN BLOCK 
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enlarged view of deformed configuration 



At limit point 



After limit point 



Element boundary overlap at the cusp 
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NEW FRICTION MODELS 


It is now the generally held view that the classical Mohr-Coloumb law of friction 
is inadequate on both physical and mathematical grounds when modelling friction 
effects in real materials. For this reason, a preliminary study of several new non- 
classical friction laws have been undertaken. These friction laws include nonlocal 
effects to approximate deformed asperities, adhesion effects, and elastic junction 
effects, and ultimately will include fracture and damage models of the junctions on 
the contact surface. The general form of these contact laws is given below, where 
CTy is the tangential frictional stress, V is the coefficient of friction, which 
may depend upon the number of loading cycles and the stress state on the contact sur- 
face, gradients of deformation, etc., S p is a smoothing operator with p represent- 
ing a characteristic dimension of deformed asperities on the contact surface, Op 
is the normal contact pressure, and <j) £ represents a compliance function which models 
the elasticity and elastoplasticity of interface junctions. This compliance reaction 
is given as an anti- symmetric function of the value of the tangential velocity vec- 
tor Uy on the contact surface. An algorithm has been developed for implementing this 
new friction law. Some preliminary results are indicated in the final two figures 
of this paper. (See also refs. 14 through 16.) j 



Oj(u) = v(Vu , t ) S (c N (u)) 4> e ( | | Uy | | ) 


! ] 


CYCLIC/DAMAGE NONLOCAL JUNCTION ELASTICITY & PLASTICITY 
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NEW FRICTION MODEL APPLICATIONS 

An example of the indentation of a rigid cylindrical stamp into an elastic 
Sl f b if 4 co nsidered . A non-classical friction law is used, with p = -0.1, 

e » E 1000, y = 0.6. Nine-node bi-quadratic elements were used. Deformed 

shape and stresses on the contact surface are shown. 
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NEW FRICTION MODEL APPLICATIONS (concluded) 
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